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ABSTRACT 

Relativistic outflows (jets) of matter are commonly observed from systems 
containing black holes. The strongest outflows occur in the radio-loud systems, 
in which the accretion disk is likely to have an advection-dominated structure. In 
these systems, it is clear that the binding energy of the accreting gas is emitted 
primarily in the form of particles rather than radiation. However, no compre- 
hensive model for the disk structure and the associated outflows has yet been 
produced. In particular, none of the existing models establishes a direct physical 
connection between the presence of the outflows and the action of a microphys- 
ical acceleration mechanism operating in the disk. In this paper we explore 
the possibihty that the relativistic protons powering the jet are accelerated at a 
standing, centrifugally-supported shock in the underlying accretion disk via the 
first-order Fermi mechanism. The theoretical analysis employed here parallels 
the early studies of cosmic-ray acceleration in supernova shock waves, and the 
particle acceleration and disk structure are treated in a coupled, self-consistent 
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manner based on a rigorous mathematical approach. We find that first-order 
Fermi acceleration at standing shocks in advection-dominated disks proves to be 
a very efficient means for accelerating the jet particles. Using physical parame- 
ters appropriate for M87 and Sgr A* , we verify that the jet kinetic luminosities 
computed using our model agree with estimates based on observations of the 
sources. 

Subject headings: accretion, accretion disks — hydrodynamics — black hole 
physics — galaxies: jets 

1. INTRODUCTION 

A large body of observational evidence has estabhshed that extragalactic relativistic jets 
are commonly associated with radio- loud active galactic nuclei (AGNs), which may contain 
hot, advection-dominated accretion disks. However, the precise nature of the mechanism 
responsible for transferring the gravitational potential energy from the infalling matter to 
the small population of nonthermal particles that escape to form the jet is not yet clear 
(see, e.g., Livio 1999). The most promising jet acceleration scenarios proposed so far arc the 
Blandford-Znajek mechanism (Blandford & Znajck 1977) and the electromagnetic cocoon 
model (Lovelace 1976; Blandford & Payne 1982), both of which involve the extraction of 
energy from the rotation of the black hole in order to power the outflow. While conceptually 
attractive, one finds that the complex physics involved in these models tends to obscure 
the nature of the fundamental microphysical processes. In particular, the introduction of 
the relativistic particles that escape to form the jet is usually made in an ad hoc manner 
without any reference to the dynamics of the associated accretion disk, although recent 
magnetohydrodynamical simulations carried out by De Villiers et al. (2005) and McKinney & 
Gammie (2004) have achieved a higher level of self-consistency. Given the relative complexity 
of the electromagnetic models, it is natural to ask whether the outflows can be explained in 
terms of well-understood microphysical processes operating in the hot, tenuous disk, such as 
the possible acceleration of the jet particles at a standing accretion shock. 

It has been known for some time that inviscid accretion disks can display both shocked 
and shock- free (i.e., smooth) solutions depending on the values of the energy and angular 
momentum per unit mass in the gas supplied at a large radius (e.g., Chakrabarti 1989a; 
Chakrabarti & Molteni 1993; Kafatos & Yang 1994; Lu & Yuan 1997; Das, Chattopadhyay, 
& Chakrabarti 2001). Shocks can also exist in viscous disks if the viscosity is relatively low 
(Chakrabarti 1996; Lu, Gu, & Yuan 1999), although smooth solutions are always possible for 
the same set of upstream parameters (Narayan, Kato, & Honma 1997; Chen, Abramowicz, 
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& Lasota 1997). Hawley, Smarr, & Wilson (1984a, 1984b) have shown through general rela- 
tivistic simulations that if the gas is falling with some rotation, then the centrifugal force can 
act as a "wall," triggering the formation of a shock. Furthermore, the possibility that shock 
instabilities may generate the quasi-periodic oscillations (QPOs) observed in some sources 
containing black holes has been pointed out by Chakrabarti, Acharyya, & Molteni (2004), 
Lanzafame, Molteni, & Chakrabarti (1998), Molteni, Sponholz, & Chakrabarti (1996), and 
Chakrabarti & Molteni (1995). Nevertheless, shocks are "optional" even when they are al- 
lowed, and one is always free to construct models that avoid them. However, in general the 
shock solution possesses a higher entropy content than the shock-free solution, and there- 
fore we argue based on the second law of thermodynamics that when possible, the shocked 
solution represents the preferred mode of accretion (Becker & Kazanas 2001; Chakrabarti & 
Molteni 1993). 

Our primary objective in this paper is to explore the consequences of the presence of 
a shock in an ADAF disk for the acceleration of the nonthermal particles in the observed 
jets. The question of whether or not viscosity needs to be included in the disk model is 
difficult to answer in general. Several authors have shown that shock solutions are possible 
in viscous (e.g., Chakrabarti 1990, 1996; Lu, Gu, & Yuan 1999; Chakrabarti & Das 2004) 
as well as inviscid disks (e.g., Chakrabarti 1989a, 1989b; Abramowicz & Chakrabarti 1990; 
Yang & Kafatos 1995; Chakrabarti 1996; Das, Chattopadhyay, & Chakrabarti 2001). In 
particular, Chakrabarti (1990) and Chakrabarti & Das (2004) demonstrated that shocks can 
exist in viscous disks if the angular momentum and the viscosity are relatively low. Since 
the acceleration of particles in shocked disks has never been investigated before, in this first 
study we shall focus on inviscid flows containing isothermal shocks (e.g., Chakrabarti 1989a; 
Kafatos & Yang 1994; Lu & Yuan 1997), while deferring the treatment of viscous disks to 
future work. However, it is clearly important to address the potential connection between 
this idealized, inviscid calculation and the physical properties of real accretion disks, which 
undoubtedly have nonzero viscosity. We argue that the results presented here should be 
qualitatively similar to those obtained in a viscous disk provided a shock is present, in which 
case efficient first-order Fermi acceleration is expected to occur. While the possible existence 
of standing shocks in viscous disks is a controversial issue at the present time, we believe 
that the work of Chakrabarti (1990, 1996), Lu, Gu, & Yuan (1999), and Chakrabarti &: Das 
(2004) provides sufficient support for the possibility to motivate the present investigation. 

Although the effect of a standing shock in heating the gas in the post-shock region 
has been examined by a number of previous authors for both viscid (Chakrabarti & Das 

2004; Lu, Gu, & Yuan 1999; Chakrabarti 1990) and inviscid (e.g., Lu & Yuan 1997, 1998; 
Yang & Kafatos 1995; Abramowicz & Chakrabarti 1990) disks, the implications of the shock 
for the acceleration of nonthermal particles in the disk have not been considered in detail 
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before. However, a great deal of attention has been focused on particle acceleration in the 
vicinity of supernova-driven shock waves as a possible explanation for the observed cosmic-ray 
energy spectrum (Blandford & Ostriker 1978; Jones & Ellison 1991). In the present paper 
we consider the analogous process occurring in hot, advection-dominated accretion flows 
(ADAFs) around black holes. These disks are ideal sites for flrst-order Fermi acceleration at 
shocks because the plasma is coUisionless and therefore a small fraction of the particles can 
gain a great deal of energy by repeatedly crossing the shock. Shock acceleration in the disk 
therefore provides an intriguing possible explanation for the powerful outflows of relativistic 
particles observed in many radio- loud systems (Le & Becker 2004). 

The dynamical model for the disk/shock/outflow employed here was discussed by Le & 
Becker (2004), who demonstrated that the predicted kinetic power in the jets agrees with the 
observational estimates for M87 and Sgr A* . Here we present a more detailed development 
of the dynamical model, including a careful examination of the implications of the shock 
acceleration process for the evolution of the relativistic particle distribution in the disk and 
the jet. The number and energy densities of the relativistic particles are determined along 
with the hydrodynamical structure of the disk in a self-consistent manner by solving the 
fluid dynamical conservation equations and the transport equation simultaneously using a 
rigorous mathematical approach. In this sense, the model presented here represents a new 
type of synthesis between studies of accretion dynamics and particle transport. 

The remainder of the paper is organized as follows. In § 2 we discuss the ADAF 
model assumptions and the possibility of shock acceleration in ADAF disks, and the general 
structure of the disk/shock model is examined in § 3. The isothermal shock jump conditions 
and the asymptotic variations of the physical parameters at both large and small radii 
are discussed in § 4. In § 5 we analyze the steady-state transport equation governing the 
distribution of the relativistic particles in the disk and the jet. Solutions for the number 
and energy density distributions of the relativistic particles are obtained in § 6, and detailed 
applications to the disks/outflows in M87 and Sgr A* are presented in § 7. The astrophysical 
implications of our results are discussed in § 8. 

2. MODEL BACKGROUND 

Accretion onto a black hole involves differentially-rotating flows in which the viscosity 
plays an essential role in transporting angular momentum outward, thereby allowing the 
accreting gas to spiral in toward the central mass (Pringle 1981). In the ADAF model, it is 
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assumed that the mass accretion rate is much smaller than the Eddington rate, 

Me = c-^ Le = 2.2 x lO"'^ (^^^ Mq yr-^ , (1) 

where the efficiency parameter (3 is of order ~ 10%, and the Eddington luminosity is defined 
by Le = 4:71 GMiTLpc/ (J^ for pure, fully-ionized hydrogen, with o"^, M, mp, and c denoting 
the Thomson cross section, the black-hole mass, the proton mass, and the speed of light, 
respectively. Due to the sub-Eddington accretion rates in these systems, the plasma is rather 
tenuous, and this strongly inhibits the efficiency of two-body radiative processes such as free- 
free emission. The gas is therefore unable to cool effectively within an accretion time, and 
consequently the gravitational potential energy dissipated by viscosity is stored in the gas as 
thermal energy instead of being radiated away (e.g., Narayan, Kato, & Honma 1997). The 
low density also reduces the level of Coulomb coupling between the ions and the electrons, 
resulting in a two-temperature configuration with the ion temperature (Tj ~ lO^^K) close 
to the virial value, and a much lower electron temperature (Tg ~ lO^K). In this scenario, 
most of the energy is advected across the horizon into the black hole, and the resulting X- 
ray luminosity is far below the Eddington value (Becker & Le 2003; Becker & Subramanian 
2005). 

When the ion temperature is close to the virial temperature, as in ADAFs, the disk is 
gravitationally unbound (e.g., Narayan, Kato, & Honma 1997; Blandford & Begelman 1999; 
Becker, Subramanian, & Kazanas 2001). It follows that the original ADAF model was not 
entirely self-consistent since it neglected outflows. This motivated Blandford & Begelman 
(1999) to propose the self-similar advection-dominated inflow-outflow solution (ADIOS) to 
address the question of self-consistency by including the possibility of powerful winds that 
carry away mass, energy, and angular momentum. In this Newtonian, nonrelativistic model, 
the dynamical solutions arc not applicable near the event horizon, and therefore the ADIOS 
approach cannot be used to obtain a global understanding of the disk structure. This led 
Becker, Subramanian, & Kazanas (2001) to modify the ADIOS scenario to include general 
relativistic effects by replacing the Newtonian potential with the pseudo-Newtonian form 
(Paczyhski & Wiita 1980) 

/ X -GM 

m - — . (2) 



s 

^2 



where = 2GM/c is the Schwarzschild radius for a black hole of mass M. This modified 
model is known as the self-similar relativistic advection-dominated infiow-outfiow solution 
(RADIOS). Despite the success of the self-similar RADIOS model in describing the general 
features of the disk/outfiow structure, it does not provide a comprehensive picture since no 
explicit microphysical acceleration mechanism is included. It is therefore natural to explore 
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possible extensions to the ADAF scenario that incorporate a concrete acceleration mechanism 
capable of powering the outflows. 

The idea of shock acceleration in the environment of AGNs was first suggested by 
Blandford & Ostriker (1978). Subsequently, Protheroe & Kazanas (1983) and Kazanas 
& EUison (1986) investigated shocks in spherically-symmetric accretion flows as a possible 
explanation for the energetic radiation emitted by many AGNs. However, in these papers the 
acceleration of the particles was studied without the benefit of a detailed transport equation, 
and the assumption of spherical symmetry precludes the treatment of acceleration in disks. 
The state of the theory was advanced by Webb & Bogdan (1987) and Spruit (1987), who 
employed a transport equation to solve for the distribution of energetic particles in a spherical 
accretion flow characterized by a self-similar velocity profile terminating at a standing shock. 
While more quantitative in nature than the earlier models, these solutions are not applicable 
to disks since the geometry is spherical and the velocity distribution is inappropriate. Hence 
none of these previous models can be used to develop a single, global, self-consistent picture 
for the acceleration of relativistic particles in an accretion disk containing a shock. 

The success of the diffusive (first-order Fermi) shock acceleration model in the cosmic- 
ray context suggests that the same mechanism may be responsible for powering the outfiows 
commonly observed in radio-loud systems containing black holes. As a preliminary step in 
evaluating the potential relevance of shock acceleration as a possible explanation for the 
observed outflows, we need to consider the basic physical properties of the hot plasma in 
ADAF disks. One of the critical issues for determining the efficiency of shock acceleration in 
accretion disks is the role of particle-particle collisions in thermalizing the high-energy ions. 
The mean free path for ion-ion collisions is given in cgs units by (Subramanian, Becker, & 
Kafatos 1996) 

where iVj and Tj denote the thermal ion number density and temperature, respectively, and 
In A is the Coulomb logarithm. In ADAF disks, Xa greatly exceeds the vertical thickness 
of the disk, and therefore the shock and the flow in general are collisionless. However, the 
mean free path Amag for collisions between ions and magnetohydrodynamical (MHD) waves 
is much shorter than for the thermal particles, and it is much longer than Xu for the 
relativistic particles (EUison & Eichler 1984; Subramanian, Becker, & Kafatos 1996). The 
increase in Amag with increasing particle energy refiects the fact that the high-energy particles 
will interact only with the highest-energy MHD waves. The low-energy background particles 
therefore tend to thermalize the energy they gain in crossing the shock due to collisions 
with magnetic waves. Conversely, the relativistic particles are able to diffuse back and forth 
across the shock many times, gaining a great deal of energy while avoiding thermalization 
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Fig. 1. — Schematic representation of our disk/shock/outflow model. The filled circles in 
the disk represent the test particles, and the unfilled circles represent the background gas. 
The test particles are injected at the shock location. 

due to the longer mean free path. 

The probability of multiple shock crossings decreases exponentially with the number of 
crossings, and the mean energy of the particles increases exponentially with the number of 
crossings. This combination of factors naturally gives rise to a power-law energy distribution, 
which is a general characteristic of Fermi processes (Fermi 1954). Two effects limit the 
maximum energy that can be achieved by the particles. First, at very high energies the 
particles will tend to lose energy to the waves due to recoil. Second, the mean free path Amag 
will eventually exceed the thickness of the disk as the particle energy is increased, resulting 
in escape from the disk without further acceleration. 



3. TRANSONIC FLOW STRUCTURE 

As discussed in § 1, various authors have established that shocks can exist in both viscid 
and inviscid disks. In this first study of particle acceleration in shocked disks, we shall focus 
on the inviscid case because it is the most straightforward to analyze from a mathematical 
viewpoint, and also because it serves to illustrate the basic physical principles involved. 
Moreover, we expect that the results obtained in the viscous case will be qualitatively similar 
to those presented here since efficient Fermi acceleration will occur whether or not viscosity 
is present, provided the fiow contains a shock. The equations governing the disk structure 
can yield solutions that include three possible types of standing shocks, namely (i) Rankine- 
Hugoniot shocks, where the effective cooling processes are so inefficient that no energy is 
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lost from the surface of the disk, (ii) isentropic shocks, where the entropy generated at the 
shock is comparable to the amount radiated away, and (ill) isothermal shocks, where the 
cooling processes are so efficient that the post-shock sound speed and disk thickness remain 
the same as the pre-shock values. In the isothermal case, the shock must radiate away both 
energy and entropy through the upper and lower surfaces of the disk (e.g., Chakrabarti 
1989a, 1989b; Abramowitz & Chakrabarti 1990). This renders the isothermal shock model 
particularly useful from the point of view of modeling outflows, since the energy lost from the 
shock can be identified with that powering the jet. On the other hand, Rankine-Hugoniot 
shocks cannot be used if we are interested in any kind of escape. The isentropic shock is an 
intermediate case. In this paper, we shall focus exclusively on the isothermal shock model 
since this case provides the strongest potential connection with the observed outfiows. 

The model considered here is depicted schematically in Figure 1. In this scenario, the gas 
is accelerated gravitationally toward the central mass, and experiences a shock transition due 
to an obstruction near the event horizon. The obstruction is provided by the "centrifugal 
barrier," which is located between the inner and outer sonic points. Particles from the 
high-energy tail of the background Maxwellian distribution are accelerated at the shock 
discontinuity via the first-order Fermi mechanism, resulting in the formation of a nonthermal, 
relativistic particle distribution in the disk. The spatial transport of the energetic particles 
within the disk is a stochastic process based on a three-dimensional random walk through 
the accreting background gas. Consequently, some of the accelerated particles diffuse to 
the disk surface and become unbound, escaping through the upper and lower edges of the 
cylindrical shock to form the outflow, while others diffuse outward radially through the disk 
or advect across the event horizon into the black hole. 

In order to analyze the connection between the disk/shock model and the transport/acceleration 
of the relativistic particles, we consider the set of physical conservation equations employed 
by Chakrabarti (1989a) and Abramowicz & Chakrabarti (1990), who investigated the struc- 
ture of a one-dimensional, steady-state, axisymmetric, inviscid accretion flow based on the 
vertically-averaged conservation equations. The effects of general relativity are incorporated 
in an approximate manner by utilizing the pseudo-Newtonian form for the gravitational 
potential per unit mass given by equation (2). The use of such a potential allows one to 
investigate the complicated physical processes taking place in the accretion disk within the 
context of a semi-classical framework while maintaining good agreement with fully relativis- 
tic calculations (see, e.g., Narayan, Kato, & Honma 1997; Becker & Subramanian 2005). 
The pseudo-Newtonian potential correctly reproduces the radius of the event horizon, the 
marginally bound orbit, and the marginally stable orbit (Paczynski & Wiita 1980). Fur- 
thermore, the dynamics of freely-falling particles near the event horizon computed using this 
potential agree perfectly with the results obtained using the Schwarzschild metric, although 



-9- 



time dilation is not included (Becker & Le 2003). 



3.1. Transport Rates 

Becker & Le (2003) and Becker & Subramanian (2005) demonstrated that three integrals 
of the flow are conserved in viscous ADAF disks, namely, the mass transport rate 

M^AnrHpv, (4) 

the angular momentum transport rate 

j = Mr^ n-g , (5) 

and the energy transport rate 

E = -gn + M(^lvl + \v^ + ^ + ^^ , (6) 

where p is the mass density, v is the radial velocity (defined to be positive for inflow), Q is 
the angular velocity, Q is the torque, H is the disk half-thickness, Vtj, — rQ, is the azimuthal 
velocity, U is the internal energy density, and P — (7 — 1)C/ is the pressure. Each of the 
various quantities represents a vertical average over the disk structure. We also assume that 
the ratio of specific heats, 7, maintains a constant value throughout the flow. Note that all 
of the transport rates M, J, and E are defined to be positive for inflow. 

The torque Q is related to the gradient of the angular velocity ft via the usual expression 
(e.g.. Prank et al. 2002) 

g^-A7rr'Hpu^ , (7) 

dr 

where is the kinematic viscosity. The disk half-thickness H is given by the standard 
hydrostatic prescription 

H{r) = ^ , (8) 
where a represents the adiabatic sound speed, defined by 



a{r) = , (9) 

and Qk denotes the Keplerian angular velocity of matter in a circular orbit at radius r in 
the pseudo- Newtonian potential (eq. [2]), defined by 

^,1.^ = 1^. (10) 
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The quantities M and J are constant throughout the flow, and therefore they represent 
the rates at which mass and angular momentum, respectively, enter the black hole. The 
energy transport rate E generally remains constant, although it will jump at the location 
of an isothermal shock if one is present in the disk. We can eliminate the torque Q be- 
tween equations (5) and (6) and combine the result with equation (9) to express the energy 
transport per unit mass as 



E 1 2 1 

— ^ - V h — + - 



e=-^^^-v'--- + ^ + - + $, (11) 



where i{r) = f2(r) and £o = J/M denote the specific angular momentum at radius r and 
the (constant) angular momentum transport per unit mass, respectively. Note that equation 
(11) is valid for both viscid and inviscid fiows. 



3.2. Inviscid Flow Equations 

In the present application, viscosity is neglected, and therefore Q = and the specific 
angular momentum is given by 

i(r) = 4 = constant (12) 

throughout the disk. It follows that the fiow is purely adiabatic, except at the location of a 
possible isothermal shock (Becker & Le 2003). In the inviscid case, equation (11) reduces to 

The resulting disk/shock model depends on three free parameters, namely the energy trans- 
port per unit mass e, the specific heats ratio 7, and the specific angular momentum L The 
value of e will jump at the location of an isothermal shock if one exists in the disk, but 
the value of £ remains constant throughout the flow. This implies that the spcciflc angular 
momentum of the particles escaping through the upper and lower surfaces of the cylindrical 
shock must be equal to the average value of the speciflc angular momentum for the particles 
remaining in the disk, and therefore the outflow exerts no torque on the disk (Becker, Sub- 
ramanian, & Kazanas 2001). Since the flow is purely adiabatic in the absence of viscosity, 
the pressure and density variations are coupled according to 

P = Dop\ (14) 



where Dq is a parameter related to the speciflc entropy that remains constant except at the 
location of the isothermal shock if one is present. 
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By combining equations (4), (8), (9), (10), and (14), we find that the quantity 

K = r3/2 {r-r^)v a^^+i)/^^-^ (15) 

is conserved throughout an adiabatic disk, except at the location of an isothermal shock. 
Following Becker & Le (2003), we refer to K as the "entropy parameter," and we note that 
the entropy per particle S is related to K via 

S^k\nK + co, (16) 

where k is the Boltzmann constant and Cq is a constant that depends on the composition of 
the gas but is independent of its state. 



3.3. Critical Conditions 



By combining equations (13) and (15), one can solve for the flow velocity v as a function 
of r using a simple root-finding procedure. However, in order to understand the implications 
of the transonic (critical) nature of the accretion flow, we must also analyze the properties of 
the "wind equation," which is the first-order differential equation governing the fiow velocity 
V. By differentiating equation (13) with respect to r, we obtain the steady-state radial 
momentum equation 

dv ^ 2a \ da , , 

(17) 



dr dr \7 — 1 J dr 

The derivative of the sound speed appearing on the right-hand side of this expression can 
be evaluated by using equations (10) and (15) to write 



1 da 
a dr 



7-1 
7 + 1 



1 dn 



K 



Qk dr 



Idv 
V dr 



(18) 



We can now construct the wind equation by combining equations (10), (17), and (18), which 
yields 

Idv _N 
V dr D ' 

where the numerator and denominator functions and D are given by 



(19) 



N 



GM 



{r-r^ 



- ^ + 



7+1 



3rg — 5r 
rir — r„) 



D 



2a^ 
7 + 1 



— V 



The simultaneous vanishing of N and D yields the critical conditions 



GM 



e 

- ^ + 



{rc-r^f rl 7 + I 



3rg — 5rc 



0, 



(20) 



(21) 
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and 



(22) 



7 + 1 

where Vc and Gc denote the values of the velocity and the sound speed at the critical radius, 



3.4. Critical Point Analysis 



Equations (21) and (22) can be solved simultaneously to express Vc and Oc as explicit 
functions of the critical radius Tc, which yields 



< = 2 



and 



(7 + 1) 



(5rc-3rg)(rc-rs)r; 



(5rc - 3rg)(rc - rg)r2 



(23) 



(24) 



The corresponding value of the entropy parameter K at the critical point is given by (see 
eq. [15]) 

K, = rf 2 (r, -r,)v, a^-^+^^Z^^"^) . (25) 

By using equations (23) and (24) to substitute for v and a in equation (13), we can express 
the energy transport parameter e in terms of Tc, i, and 7, obtaining 



' 2 r2 



GM 



27 
7-1 



GMr^-(Hr,^-rJ' 
'^c('^c-rs)(5rc-3rs) 



This expression can be rewritten as a quartic equation for Tc of the form 



(26) 



(27) 



where 



12 e + 



8 



1/5-7 



2 V7- 1 



£2-6, 



^f = 5e, C)=16e-3 + 

Q 
7^ 



7-1 



(28) 



7-1. 
27 + 6\ ,2 



7-1 
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and we have utilized natural gravitational units with GM — c — 1 and — 2. These equa- 
tions agree with the corresponding results derived by Das, Chattopadhyay, & Chakrabarti 
(2001). The four solutions for Vc in terms of the three fundamental parameters e, i, and 7 can 
be obtained analytically using the standard formulas for quartic equations (e.g., Abramowitz 
& Stegun 1970). We refer to the roots using the notation Vd, rc2, Tcs, and Tca in order of 
decreasing radius. 

The critical radius always lies inside the event horizon and is therefore not physically 
relevant, but the other three are located outside the horizon. The type of each critical 
point is determined by computing the two possible values for the derivative dv/dr at the 
corresponding location using L'Hopitars rule and then checking to see whether they are real 
or complex. We find that both values are complex at rc2, and therefore this is an 0-type 
critical point, which does not yield a physically acceptable solution. The remaining roots Tci 
and each possess real derivatives, and are therefore physically acceptable sonic points, 
although the types of accretion flows that can pass through them are different. Specifically, 
rc3 is an X-type critical point, and therefore a smooth, global, shock-free solution always 
exists in which the flow is transonic at rcs and then remains supersonic all the way to 
the event horizon. On the other hand, r^i is an «-type critical point, and therefore any 
accretion flow originating at a large distance that passes through this point must display 
a shock transition below Td (Abramowicz & Chakrabarti 1990). After crossing the shock, 
the subsonic gas must pass through another (a-type) critical point before it enters the black 
hole since the flow has to be supersonic at the event horizon (Weinberg 1972). 



3.5. Shock-Free Solutions 

Even when a shock can exist in the flow, it is always possible to construct a globally 
smooth (shock-free) solution using the same set of parameters. Smooth flows must pass 
through the inner critical point located at radius Tcs, which is calculated using the quartic 
equation (27) for given values of e, I, and 7. The corresponding values for the critical 
velocity, Vcs, the critical sound speed, Ocs, and the critical entropy, ifcs, are computed using 
equations (23), (24), and (25), respectively. Since the flows treated here are inviscid, they 
have a conserved value for the entropy parameter K (eq. [15]) unless a shock is present. 
Hence in a smooth, shock-free flow, the value of K is everywhere equal to the critical value 
Kf^. The structure of the velocity profile in a shock-free disk can therefore be determined 
using a simple root-finding procedure as follows. By eliminating the sound speed a between 



equations (13) and (15), we obtain the equivalent expression 




1 (7-l)/{7+l) 



(29) 



where we have set K = K^^. In general, at any radius r, equation (29) yields one subsonic root 
and one supersonic root for the velocity. The subsonic solution is chosen for r > r^s, and the 
supersonic solution is selected for r < r^a. Once the velocity profile v{r) has been computed, 
wc can obtain the corresponding sound speed distribution a(r) by utilizing equation (15) 
with K = Kc3. Note that the velocity and sound speed solutions can also be calculated by 
integrating the wind equation (19) numerically, and the results obtained using this approach 
agree with the root-finding method. 



Our primary goal in this paper is to analyze the acceleration of relativistic particles due 
to the presence of a standing, isothermal shock in an accretion disk. Hence we are interested 
in fiows that pass smoothly through the outer critical radius rd, and then experience a 
velocity discontinuity at the shock location, which we refer to as r*. In order to form self- 
consistent global models, we will first need to understand how the structure of the disk 
responds to the presence of a shock. This requires analysis of the shock jump conditions, 
which are based on the standard fiuid dynamical conservation equations. Since shocks are 
always optional even when they can occur, we will compare our results for the relativistic 
particle acceleration with those obtained when there is no shock and the fiow is globally 
smooth. 

The values of the energy transport parameter e on the upstream and downstream sides 
of the isothermal shock at r = are denoted by e_ and e_|_, respectively. Note that e_ > e_|_ 
as a consequence of the loss of energy through the upper and lower surfaces of the disk at the 
shock location. It is important to emphasize that the drop in e at the shock has the effect 
of altering the transonic structure of the fiow in the post-shock region. Hence, although 
the post-shock fiow must pass through another critical point and become supersonic before 
crossing the event horizon, the new (inner) critical point is not equal to any of the four roots 
computed using the upstream energy transport parameter e_. Instead, the new inner critical 
radius, which we refer to as fcs, must be computed using the downstream value of the energy 
transport parameter, e+. We point out that the total energy inflow rate across the horizon, 
including the rest-mass contribution, must be positive since no energy can escape from the 
black hole, and therefore we require that -|- e_|_ > 0. 



4. 



ISOTHERMAL SHOCK MODEL 
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4.1. Isothermal Shock Jump Conditions 

We shall assume that the escape of the relativistic particles from the disk results in a 
negligible amount of mass loss because the Lorentz factor of the escaping particles is much 
greater than unity (see Table 2). This is confirmed ex post facto by comparing the rate of 
mass loss, Mioss, with the accretion rate M. We find that for the models analyzed here, 
M\oss/M < 10~^, and therefore our assumption of neghgible mass loss is justified. Hence 
the accretion rate M is conserved as the gas crosses the shock, which is represented by the 
condition 

AM= limM{n-e)- M{n + e) = 0, (30) 

where the symbol "A" will be used to denote the difference between post- and pre-shock 
quantities. The specific angular momentum i = J / M is also conserved throughout the flow, 
and therefore we find that 

A j = . (31) 

Furthermore, the radial momentum transport rate, defined by 

i = A7rrH(P + pv'^) , (32) 

must remain constant across the shock, and consequently 

A/ = 0. (33) 

Based on equations (13) and (31), we find that the jump condition for the energy transport 
rate E is given by 

A^ = M Qa^;2 + -1^ Aa^^ . (34) 

Equations (4), (8), and (13) can be combined with equations (30), (33), and (34) to 
obtain 

p+f+a+ = /0_f^a_ , (35) 
a+P+ — tt-P- = a^p^v'i — a+p+f ^ , (36) 
vl -v'i ai - a'i , . 

where the subscripts "-" and "+" refer to quantities measured just upstream and just down- 
stream from the shock, respectively. In the case of an isothermal shock, a+ = a_, and 
therefore the shock jump conditions reduce to 

P+V+ = p-v- , (38) 
P+-P- = p-v^_-p+vl, (39) 

e+-e_ = (40) 
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Fig. 2. — Plot of the (e_,^) parameter space for an ADAF disk with 7 = 1.5. Only smooth 
flows exist in region I, and both shocked and smooth solutions are possible in regions II and 
III. When i > ^max (region IV), no steady-state dynamical solutions can be developed. The 
parameters corresponding to models 2 and 5 are indicated. 

Combining equations (38) and (39) and substituting for the pressure P using equation (9) 
yields the velocity jump condition 

^^^-^MZ' < 1, M- = — , (41) 

where Al_ is the incident Mach number of the shock. The corresponding result for the shock 
compression ratio i?* is 

R, = P±^^Ml > 1. (42) 
P- 

Hence the gas density increases across the shock as expected. Based on equations (15) and 
(41) and the fact that a+ = a_ in an isothermal shock, we find that the jump condition for 
the entropy parameter K is given by 

^^r'M-J < 1, (43) 

which indicates that entropy is lost from the disk at the shock location due to the escape of 
the particles that form the outflow (jet). 
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We can also make use of equation (41) to rewrite the jump condition for the energy 
transport parameter (eq. [40]) as 



Ae^e,-e_ = -^-^-lJ < 0. (44) 

The associated rate at which energy escapes from the disk at the isothermal shock location 
(the "shock luminosity") is given by (see eq. [13]) 

L,hock = -A£; = -MAe ocergss-^. (45) 

Eliminating Ae between equations (44) and (45) yields the alternative result 

= '!L (l (46) 



4.2. Shock Point Analysis 

For a given value of 7, it is known that smooth, shock-free global flow solutions exist only 
within a restricted region of the (c-,^) parameter space, and isothermal shocks can occur 
only in a subset of the smooth-flow region. In order for a shock to exist in the flow, it must 
be located between two critical points, and it must also satisfy the jump conditions given 
by equations (41), (43), and (44). The procedure for determining the disk/shock structure 
is summarized below. 

The process begins with the selection of values for the fundamental parameters e_, £, 
and 7. The values of e_ and ^ are ultimately constrained by the observations of a specific 
object, as discussed in § 7. Following Narayan, Kato, & Honma (1997), we shall assume 
an approximate equipartition between the gas and magnetic pressures, and therefore we set 
7 = 1.5. The first step in the determination of the shock location is the computation of 
the outer critical point location, rd, using the quartic equation (27). The associated values 
for the critical velocity, Vci, the critical sound speed, Od, and the critical entropy, are 
then calculated using equations (23), (24), and (25), respectively. Note that since the fiow 
is adiabatic everywhere in the pre-shock region, it follows that 

K_ = . (47) 

The profiles of the velocity v{r) and the sound speed a(r) in the pre-shock region can therefore 
be calculated using a root-finding procedure based on equation (29), or, alternatively, by 
integrating numerically the wind equation (19). The next step is the selection of an initial 



guess for the shock radius, r*, and the calculation of the associated shock Mach number 
M.- = u_/a_ using the pre-shock dynamical solutions for v{r) and a(r). Based on the value 
of we can compute the jump in the entropy parameter K using equation (43), and 

consequently we find that the entropy in the downstream region is given by 

where we have also used equation (47) . 

In order to determine whether the initial guess for r^, is self-consistent, we employ a 
second, independent procedure for calculating the entropy in the downstream region based 
on the critical nature of the flow. In this approach, the downstream energy parameter e+ is 
calculated using the jump condition given by equation (44), which yields 

We utilize this value to compute the downstream critical point radius based on the quartic 
equation (27). The associated values for the critical velocity, Vcs, the critical sound speed, 
Ocs, and the critical entropy, K^s, are then calculated using equations (23), (24), and (25), 
respectively. The final step is to compare the value of Kc3 with that obtained for using 
equation (48). If these two quantities are equal, then the shock radius r* is correct and the 
disk/shock model is therefore dynamically self-consistent. Otherwise, the value for must 
be iterated and the search continued. Roots for can be found only in certain regions of 
the (e_,£, 7) parameter space, as discussed below. 

By combining the analysis of the shock location discussed above with the critical con- 
ditions developed in § 3, we are able to compute the structure of shocked and shock-free 
(smooth) disk solutions for a given set of parameters e_, i, and 7. The resulting topology 
of the parameter space is depicted in Figure 2 for the case with 7 = 1.5, which is the main 
focus of this paper. Within region I, only smooth flows are possible, and in regions II and III 
both smooth and shocked solutions are available. No global flow solutions (either smooth or 
shocked) exist in region IV, with i > £max- Inside region II, one root for the shock radius r* 
can be found, and in region III two shock solutions are available, although only one actual 
shock can occur in a given flow. It is unclear which of the two roots for in region III 
is preferred since the stability properties of the shocks are not completely understood (e.g., 
Chakrabarti 1989a, 1989b; Abramowicz & Chakrabarti 1990). However, it is worth noting 
that the inner shock is always the stronger of the two because the Mach number diverges 
as the gas approaches the horizon. The larger compression ratio associated with the inner 
location leads to more efficient particle acceleration and enhanced entropy generation, and 



(48) 
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therefore we expect that the inner solution is preferred in nature. Based on this argument, 
we will focus on the inner shock location in our subsequent analysis. 

Before we proceed to examine the transport equation for the relativistic particles, it is 
important to analyze the asymptotic solutions obtained for the dynamical variables near the 
event horizon and also at large radii. This is a crucial issue since the nature of the global 
solutions to the transport equation depends sensitively on the boundary conditions imposed 
at large and small radii. The asymptotic solutions for the dynamical variables near the event 
horizon were fully discussed by Becker & Le (2003). We shall briefly review their results 
and then perform a similar analysis in order to determine the asymptotic properties of the 
solutions as r — oo. 



4.3. Asymptotic Behavior Near the Horizon 

Becker & Le (2003) demonstrated that the variation of the global solutions in a viscous 
ADAF disk becomes purely adiabatic close to the event horizon, and therefore the asymp- 
totic solutions that they obtain can be directly applied to our inviscid model. Using their 
equations (47) and (51), we find that the asymptotic variations of the radial velocity v and 
the sound speed a near the horizon are given by 

v^r) oc (r - r^)-' , a^r) oc (r - r^Y'-^y('+^') , r ^ . (50) 

The divergence of as r — > Tg implies that it cannot represent the standard velocity in 
the region near the horizon. However, our dynamical model is consistent with relativity 
if we interpret v as the radial component of the four-velocity in this region (Becker & Le 
2003; Becker & Subramanian 2005). By combining these relations with equations (4), (8), 
and (10), we find that the corresponding results for the asymptotic variations of the disk 
half-thickness H and the density p become 

H{r) oc (r - r^)(.y+mw) ^ p(^) _ rJ-i/(7+i) ^ r ^ r, . (51) 

4.4. Asymptotic Behavior at Infinity 

We can use the energy transport equation (13) and the entropy equation (15) to obtain 
the asymptotic solutions for v and a at infinity as follows. In the limit r — > oo, the two 
dominant terms in equation (13) are e and 0^/(7 — 1), and therefore we find that 



a^(r) — > (7 — 1) e , r — > 00 



(52) 
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Recalling that K is constant in the adiabatic upstream flow, we can combine equations (15) 
and (52) to conclude that the asymptotic variation of the inflow velocity is given by 

V oc r~^/^ , r — >• oo . (53) 

Finally, based on equations (8), (10), and (52), we flnd that the disk half-height varies as 

H oc r^/^ , r ^ oo . (54) 

We can also combine equations (9), (14), and (52) to conclude that the asymptotic behavior 
of the density is given by 

p — > const. , r — > oo . (55) 



5. STEADY-STATE PARTICLE ACCELERATION 

Our goal in this paper is to analyze the transport and acceleration of relativistic particles 
(ions) in a disk governed by the dynamical model developed in § 3-4. For flxed values of 
the theory parameters e_, i, and 7, we will study the transport of particles in disks with 
and without shocks. The particle transport model utilized here includes advection, spatial 
diffusion, Fermi energization, and particle escape. In order to maintain consistency with 
the dynamics of the disk, we will need to equate the energy escape rate for the relativistic 
particles with the "shock luminosity" Lghock given by equation (45). Our treatment of Fermi 
energization includes both the general compression related to the overall convergence of 
the accretion flow, as well as the enhanced compression that occurs at the shock. In the 
scenario under consideration here, the escape of particles from the disk occurs via vertical 
spatial diffusion in the tangled magnetic fleld, as depicted in Figure 1. To avoid unnecessary 
complexity, we will utilize a simplified model in which only the radial (r) component of the 
spatial particle transport is treated in detail. In this approach, the diffusion and escape of 
the particles in the vertical [z) direction is modeled using an escape-probability formalism. 
We will treat the relativistic ions as test particles, meaning that their contribution to the 
pressure in the flow is neglected. This assumption is valid provided the pressure of the 
relativistic particles turns out to be a small fraction of the thermal pressure, as discussed in 
§ 8. The ions accelerated at the shock are energized via collisions with MHD waves advected 
along with the background (thermal) flow, and therefore the shock width is expected to be 
comparable to the magnetic coherence length, Amag- This approximation will be used to 
determine the rate at which particles escape from the disk in the vicinity of the shock. 
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5.1. Transport Equation 

The Green's function, f^{Eo, E, r^,,r), represents the particle distribution resulting from 
the continual injection of A^o particles per second, each with energy Eo, from a source located 
at the shock radius, r = r^.. In a steady-state situation, the Green's function satisfies the 
transport equation (Becker 1992) 

^ = = -V-F-^^ (^E'^ . v/^) + _ /^^^ , (56) 

— * 

where the specific fiux F is evaluated using 

and the source and escape terms are given by 

• No5{E-Eo)5{r-r.) ■ 

^ {inE^rnH. ' ^ ' ~ ■ ^^^^ 

The quantities E, k, v, and = H{r^) represent the particle energy, the spatial diffusion 
coefficient, the vector velocity, and the disk half-thickness at the shock location, respectively, 
and the dimensionless parameter Aq determines the rate of particle escape through the 
surface of the disk at the shock location. The vector velocity v has components given by 
V = Vrf + VzZ + f (j), where v = —Vr is the positive inflow speed. 

The total number and energy densities of the relativistic particles, denoted by and 
Ur, respectively, are related to the Green's function via 

/•oo roo 

nr{r)= ATvE'^f^dE, Ur{r) = j A-kE^ f^dE , (59) 
Jo Jo 

which determine the normalization of f^. Equations (56) and (57) can be combined to obtain 
the alternative form 

^ • V/o = ^ £^ H- + V • V/^) + /source - /esc , (60) 

where the left-hand side represents the co-moving (advective) time derivative and the terms 
on the right-hand side describe flrst-order Fermi acceleration, spatial diffusion, the particle 
source, and the escape of particles from the disk at the shock location, respectively. Note 
that escape and particle injection are localized to the shock radius due to the presence of 
the 5-functions in equations (58). Our focus here is on the first-order Fermi acceleration of 
relativistic particles at a standing shock in an accretion disk, and therefore equation (60) 
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does not include second-order Fermi processes that may also occur in the flow due to MHD 
turbulence (e.g., Schhckeiser 1989a,b; Subramanian, Becker, & Kazanas 1999). 

Under the assumption of cylindrical symmetry, equations (58) and (60) can be rewritten 

as 



Id. , dvz 
r or dz 



dr dz 3 

NoS{E-Eo) S{r-r 



E—^ — rK 



{AnEoyr.H, 



dE r dr \ dr 

Aoc6{r-r.)f^ , (61) 



where the escape of particles from the disk is described by the final term on the right-hand 
side. In Appendix A, we demonstrate that the vertically integrated transport equation is 
given by (sec cq. [A9]) 



dr 3r dr dE r dr \ dr 

NoS{E-Eo)5{r-r.) , „ , . .... 

+ ^^^^^ A,cH.f,S{r-r.), (62) 

where the symbols f^, Vr, and k represent vertically averaged quantities. We establish in 
Appendix B that within the context of our one-dimensional spatial model, the dimensionless 
escape parameter appearing in equation (62) is given by (see eqs. [B8] and [BIO]) 

o \ 2 



where = («;_ -|- k,+)/2 denotes the mean value of the diffusion coefficient at the shock 
location. The condition < 1 is required for the validity of the diffusive picture we have 
employed in our model for the vertical transport. 



5.2. Number and Energy Densities 

The energy moments of the Green's function, /n(?"), are defined by 

4(r)= / inE^f^dE, 
Jo 

so that (cf. eqs. [59]) 



(64) 



nr{r) = hir) , C/,(r) = h{r) . (65) 
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By operating on equation (62) with AnE" dE and integrating by parts once, we find that 
the function 7„ satisfies the differential equation 

^ dr \ 3 J r dr ^ r dr \ dr J 

+ -AcH.U{r-r.), (66) 



which can be expressed in the flux-conservation form 

S(r — r .) 

(67) 



4- UTrrHFn) = AurH 
dr 



^-n\ dig NoEr'5{r-r.) 

""dr^ Anr.H. "^'"^^"^ 



where 47rrifF„ represents the rate of transport of the nth moment, and the flux F„ is defined 
by 

n + 1 \ ^ din 



and V — —Vr denotes the positive infiow speed. 

In order to close the system of equations and solve for the relativistic particle number 
and energy densities hir) and hir) using equation (66), we must also specify the radial 
variation of the diffusion coefficient k. The behavior of k can be constrained by considering 
the fundamental physical principles governing accretion onto a black hole. First, we note 
that near the event horizon, particles are swept into the black hole at the speed of light, 
and therefore advcction must dominate over diffusion. This condition applies to both the 
thermal (background) and the nonthermal (relativistic) particles. Second, we note that in 
the outer region (r —>■ oo), diffusion is expected to dominate over advection. Focusing on the 
flux equation for the particle number density rir, obtained by setting n = 2 in equation (68), 
we can employ scale analysis to conclude based on our physical constraints that 

lini ^ ^ ^ = , hm ^ = . (69) 

r^rg [r — r^) V{r) r-^oo K.[r) 

The precise functional form for the spatial variation of k is not completely understood in the 
accretion disk environment. In order to obtain a mathematically tractable set of equations 
with a reasonable physical behavior, we shall utilize the general form 

K{r) ^ Kov{r)rs -Ij , (70) 



where ko and a are dimensionless constants. Due to the appearance of the inflow speed 
V in equation (70), we note that k exhibits a jump at the shock. This is expected on 
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physical grounds since the MHD waves that scatter the ions are swept along with the thermal 
background flow, and therefore they should also experience a density compression at the 
shock. 

As discussed above, close to the event horizon, inward advection at the speed of hght 
must dominate over outward diffusion. Conversely, in the outer region, we expect that 
diffusion will dominate over advection as r — > oo. By combining equation (70) with the 
asymptotic velocity variations expressed by equations (50) and (53), we find that the con- 
ditions given by equations (69) are satisfied if « > 1, and in our work we set a = 2. Note 
that the escape parameter Aq is related to kq via equation (63), which can be combined with 
equation (70) to write 



where = {v^ + v+)/2 represents the mean velocity at the shock location r = r*. The 
value of the diffusion parameter kq is constrained by the inequality in equation (71). In § 7 
we demonstrate that ko can be computed for a given source based on energy conservation 
considerations. With the introduction of equations (70) and (71), we have completely defined 
all of the quantities in the transport equation, and we can now solve for the number and 
energy densities of the relativistic particles. The particle distribution Green's function 
and its applications will be discussed in a separate paper. 



Once the disk/shock dynamics have been computed based on the selected values for 
the free parameters e_, £, and 7 using the results in §§ 3 and 4, the associated solutions for 
the number and energy densities of the relativistic particles in the disk can be obtained by 
solving equation (66). In the case of the number density, n,, = I2, an exact solution can 
be derived based on the linear first-order differential equation describing the conservation of 
particle flux. However, in order to understand the variation of the energy density, Ur = I3, 
we must numerically integrate a second-order equation. The solutions obtained below are 
applied in § 7 to model the outflows observed in M87 and Sgr A* . 




(71) 



6. SOLUTIONS FOR THE ENERGY MOMENTS 
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6.1. Relativistic Particle Number Density 

The equation governing the transport of the particle number density, n^, is obtained by 
setting n = 2 in equation (67), which yields 

dNr 

—!- = Nq 5{r - r*) - ^.tit^H^Aq cn^ Sir - r\) , (72) 
or 

where the relativistic particle transport rate, Nr{r), is defined by (cf. eq. [68]) 

Nrir) = -AtttH {vUr + Ki oc s"^ , (73) 

and Nj. > Q ii the transport is in the outward direction. Since the source is located at the 
shock, there are two spatial domains of interest in our calculation of the particle transport, 
namely domain I (r > r*), and domain II (r < r*). Note that the number density nr{r) must 
be continuous at r = r* in order to avoid generating an infinite diffusive flux according to 
equation (73). Away from the shock location, r 7^ r*, and therefore equation (72) reduces 
to Nr = const., which reflects the fact that particle injection and escape are localized at the 
shock. We can therefore write 

N(r) = I ' ^ (74) 

where the constant TVi > denotes the rate at which particles are transported outward, 
radially, from the source location, and the constant Nu < represents the rate at which 
particles are transported inward towards the event horizon. 

The magnitude of the jump in the particle transport rate at the shock is obtained by 
integrating equation (72) with respect to radius in a very small region around r — r^, which 
yields 

N, - Nn ^No- A^esc , (75) 

where 

A^esc = ^nr^H^Aocn^ (76) 

represents the (positive) rate at which particles escape from the disk at the shock location 
to form the outflow (jet), and = ^^(r*). If no shock is present in the flow, then Aq — 
and therefore iVcsc=0. Note that the discontinuity in Nr at the shock produces a jump in 
the derivative drir/dr via equation (73). 

We can rewrite equation (73) for the number density in the form 

dr K ^ AnrHn ' 



-26- 



which is a linear, first-order differential equation for nr{r). Using the standard integrating 
factor technique and employing equation (70) for k yields the exact solution 



nr{r) — e 



-J{r) 



Nrjr) 
An 



r'HK 



dr' 



(78) 



where Nr{r) is given by equation (74) and the function J(r) is defined by 

.M./>.-'d/^.„-[(^-i)-'-(^-i)-' 



(79) 



According to equation (78), nr{r) is continuous at the shock/source location as required. 
Far from the black hole, diffusion dominates the particle transport, and therefore should 
vanish as r — > oo. In order to ensure this behavior, we must have 



where 



Att .L r'HK 



dr' . 



(80) 



(81) 



Furthermore, in order to avoid exponential divergence of as r — > r^ in domain II, we also 
require that 

= -7V„C„, (82) 

where 



Cn = 



1 r* e^^''') 



47r 



r'HK 



dr' 



(83) 



By combining equations (75), (76), (80), and (82), we can develop explicit expressions 
for the quantities n*, Ni, Nu, and iVesc based on the values of r* and Nq and the profiles of 
the inflow velocity v{r) and the diffusion coefficient K{r). The results obtained are 



Cf ^ + QT^ + Airr^H^Ao c 



(84) 



^11 = --FT, 
^11 



Anr^H^Aocn^ . 



These relations, along with equation (78), complete the formal solution for the relativistic 
particle number density ^^(r). The solution is valid in both shocked and shock- free disks 
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(the shock- free case is treated by setting Aq — 0). When a shock is present, the particle 
escape rate TVesc is proportional to Nq but is independent of Eq by virtue of equations (84). 

It is interesting to examine the asymptotic variation of rir near the event horizon and 
also at large distances from the black hole. Far from the hole, advection is neghgible and 
the particle transport in the disk is dominated by outward-bound diffusion. In this case we 
can use equation (77) to conclude that 

drir Ni , . 

-1 7]-: r^oo, (85 

ar AirrHK 

where we have used the fact that A^^ = Ni for r > r*. By combining equations (70) and (85) 

with the asymptotic relations given by equations (53) and (54), we find upon integration 
that ^ 

nr{r) oc - , r ^ oo . (86) 

In order to study the behavior of rir near the event horizon, we take the limit as r — >■ Tg in 
equation (78), obtaining after some algebra 

^ -4^ ' ' ^ ^^^^ 
where we have set Nr — Nu- Comparing this relation with equation (4), we find that 

nr{r) oc p(r) , r — > Tg , (88) 

where p is the density of the background (thermal) gas. Equation (88) is a natural conse- 
quence of the fact that the particle transport near the horizon is dominated by inward-bound 
advection. We can also combine equations (51) and (88) to obtain the explicit asymptotic 
form 

nr{r) oc (r - rg)-^/(^+^) , r ^ . (89) 



6.2. Relativistic Particle Energy Density 

The differential equation satisfied by the relativistic particle energy density, Ur — h, is 
obtained by setting n = 3 in equation (66), which yields 



dr 3 r dr r dr \ dr 

+ ^^^^^^^^^-AcH.UrSir-r.). (90) 
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By analogy with equations (72) and (73), we can recast this expression in the flux-conservation 
form 

dEr 



dr 



= 47rri7 



vdUr NoEo5{r-n) s(r-r) 



(91) 



where the relativistic particle energy transport rate, Er{r), is defined by 



Er{r) = -AtttH {^vUr + K (X ergs s~^ , (92) 

and Er > Q for outwardly directed transport. Note that unlike the number transport rate 
iVj., the energy transport rate E^ does not remain constant within domains I and II due to 
the appearance of the first term on the right-hand side of equation (91), which expresses the 
compressional work done on the relativistic particles by the background flow. 

Although the energy density U,. must be continuous at the shock/source location in order 
to avoid generating an infinite diffusive ffux, the derivative dUr/dr displays a discontinuity 
a,t r — r*, which is related to the jump in the energy transport rate via equation (92). By 
integrating equation (91) in a very small region around r = r*, we find that 

AEr^L^sc-NoEo , (93) 

where 

Lesc = 47rr*if*Ao c C/* oc ergs s"^ (94) 

denotes the rate of escape of energy from the disk into the outflow (jet) at the shock location, 
and [/* = Ur{r^). If no shock is present, then = and therefore Lesc = 0. We remind the 
reader that the symbol "A" refers to the difference between post- and pre-shock quantities 
(see eq. [30]). Equations (92) and (93) can be combined to show that the derivative jump is 
given by 

dUr\ iVo^O-Lesc 4 

The differential equation (90) governing the relativistic particle energy density is second- 
order in radius, and therefore we will need to establish two boundary conditions in order to 
solve for Ur{r). These can be obtained by analyzing the behavior of Ur close to the event 
horizon and at large distances from the black hole. Far from the hole, advection is negligible 
and the particle transport in the disk is dominated by outward-bound diffusion. In this 
regime, Fermi acceleration is negligible, and consequently we find that Ur oc n^. We can 
therefore use equation (86) to conclude that 

Ur{r) (X - , r — > oo . (96) 
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Close to the event horizon, the particle transport is dominated by advection, and therefore 
Ur and rir obey the standard adiabatic relation 

Ur oc , r ^ Tg . (97) 

Combining this result with equations (88) then yields 

UrOc{r- rg) -4/(3^+3) , r ^ . (98) 

The global solution for Ur{r) can now be expressed as 



where A and B are constants and the functions Qi{r) and Qii{r) satisfy the homogeneous 
differential equation (see eq. [90]) 

H./^ = -i^^irm^)^'-^(rHj-f] . (100) 



dr 3 r dr r dr \ dr 

along with the boundary conditions (see eqs. [96] and [98]) 

/r /r \ -4/(37+3) 

Qiirout) ^ [y^j , Qii(rin) = (^^ - ij , (101) 

where r^^ and rout denote the radii at which the inner and outer boundary conditions are ap- 
plied, respectively. The constants A and B arc computed by requiring that Ur be continuous 
at r = r* and that the derivative dUr/dr satisfy the jump condition given by equation (95). 
The results obtained are 



B 



b9r 

Qi 

NqEq 
ATir^H^ 



- {V_ - V+) Qu + QuK'+ h AqcQu 

6 LJi 



(102) 
(103) 



where primes denote differentiation with respect to radius. The solutions for the functions 
Qi{r) and Qii{r) are obtained by integrating numerically equation (100) subject to the 
boundary conditions given by equations (101). Once the constants A and B are computed 
using equations (102) and (103), the global solution for Ur{r) is evaluated using equation (99). 
The solution applies whether or not a shock is present in the flow. The shock-free case is 
treated by setting Aq — 0. This completes the solution procedure for the relativistic particle 
energy density. The results derived in this section are used in § 7 to model the outflows 
observed in M87 and Sgr A* . 
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7. ASTROPHYSICAL APPLICATIONS 

Our goal is to determine the properties of the integrated disk/shock/outflow model 
based on the observed values for the black hole mass M, the mass accretion rate M, and 
the jet kinetic power Ljet associated with a given source. The fundamental free parameters 
for the theoretical model are e_, £, and 7. Since we set 7 = 1.5 in order to represent an 
approximate equipartition between the gas and magnetic pressures (e.g., Narayan, Kato, & 
Honma 1997), only e_ and i remain to be determined. Here we describe how global energy 
conservation considerations can be used to solve for the various theoretical parameters in 
the model based on observations. 



7.1. Energy Conservation Conditions 

Once the values of M, M, and Ljet have been specified for a source based on observations, 
we select a value for the free parameter e_ and then compute £ by satisfying the relation 

-t'shock = L^et , (104) 

where I/ghock is the shock luminosity given by equation (46). This result ensures that the 
jump in the energy transport rate at the isothermal shock location is equal to the observed 
jet kinetic luminosity. The procedure for determining i also includes solving for the shock 
location and the critical structure using results from §§3 and 4. The velocity profile v{r) is 
computed either by numerically integrating equation (19) or by using a root-finding proce- 
dure based on equation (29), and the associated solution for the adiabatic sound speed a(r) 
is obtained using equation (13). 

After the velocity profile has been determined, we can compute the number and energy 
density distributions for the relativistic particles in the disk using equations (78) and (99), 
respectively. This requires the specification of the injection energy of the seed particles 
Eq as well as their injection rate Nq. We set the injection energy using Eq = 0.002 ergs, 
which corresponds to an injected Lorentz factor Tq = EQ/^rrLpC^) ~ 1.3, where rrip is the 
proton mass. Particles injected with energy Eq are subsequently accelerated to much higher 
energies due to repeated shock crossings. We find that the speed of the injected particles, 
vo — c{l — Fq ^)^/^, is about three to four times higher than the mean ion thermal velocity at 
the shock location, Vrms — (3A;T*/mp)^/^, where T* is the ion temperature at the shock. The 
seed particles are therefore picked up from the high-energy tail of the MaxweUian distribution 
for the thermal ions. With £^0 specified, we can compute the particle injection rate Nq using 
the energy conservation condition 

TVo ^0 = , (105) 
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Table 1: Disk structure parameters 

Model e_ £ e_|_ TcI ^£3 ^c3 jVl» R* 

1 0.001600 3.1040 -0.005676 93.177 5.478 19.624 5.798 10.369 1.094 1.795 1.28 

2 0.001527 3.1340 -0.005746 98.524 5.379 21.654 5.659 11.544 1.125 1.897 1.16 

3 0.001400 3.1765 -0.005875 109.781 5.252 24.500 5.487 13.154 1.170 2.052 1.01 

4 0.001240 3.1280 -0.008723 131.204 5.408 14.073 5.850 6.819 1.113 1.857 1.65 

5 0.001229 3.1524 -0.008749 131.874 5.329 15.583 5.723 7.672 1.146 1.970 1.49 

6 0.001200 3.1756 -0.008778 135.192 5.260 16.950 5.614 8.434 1.177 2.076 1.36 

All quantities are expressed in gravitational units {GM = c = 1) except T* which is in units of 10" K. 

which ensures that the rate at which energy is injected into the flow in the form of the 
relativistic seed particles is equal to the energy loss rate for the background gas at the 
isothermal shock location. 

In order to maintain agreement between the transport model and the observations, we 
must also require that the rate at which particle energy escapes from the disk due to vertical 
diffusion is equal to the observed jet power. This condition can be written as 

-^esc = -^jet , (106) 

where L^sc is the energy escape rate given by equation (94). The escape constant, Aq, 
appearing in the transport equation is independent of the particle energy in our model, and 
consequently the escaping particles will have exactly the same mean energy as those in the 
disk at the shock location. The mean energy of the escaping particles is therefore given by 

^esc = — , (107) 

where and f/* denote the number and energy densities of the relativistic particles at the 
shock location, respectively. Hence Ecsc is proportional to Eq but it is independent of A^o- 
We note that equations (76), (94), and (107) can be combined to show that 

-^esc — Nesc -^esc , (108) 

where TVesc is the particle escape rate (eq. [76]). By satisfying equations (104), (105), and 
(106), we ensure that energy is properly conserved in our model. Taken together, these 
relations allow us to solve for the various theoretical parameters based on observational 
values for M, M, and Ljet, as explained below. 



7.2. Model Parameters 



Our simulations of the disk structure and particle transport in M87 and Sgr A* are 
based on various published observational estimates for M, M, and Ljet- For M87, we set 
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Table 2: Transport equation parameters 



Model 


No 


Ni 


Ko 




Ao 




U* 


Neec 


-^esc 


Tesc 




(s-i) 


Nu 








(cm-3) 


(ergs cm-3) 


No 


Eq 




2 


2.75 xlO'*** 


-0.18 


0.02044 


0.427877 


0.0124 


2.01 X 10^ 


2.39 X 10^ 


0.17 


5.95 


7.92 


5 


2.51 xlO*i 


-0.15 


0.02819 


0.321414 


0.0158 


4.33 X 10^ 


4.71 X 103 


0.18 


5.45 


7.26 



All quantities are expressed in gravitational units {GM = c = 1) except as noted. 

M = 3 X 10^ Mq (e.g., Ford et al. 1994), M = 1.3 x 10"^ Mq yr'^ (e.g., Reynolds et al. 
1996), and Ljet = 5.5 x 10'^^ ergs s~^ (e.g., Reynolds et al. 1996; Bicknell & Begelman 1996; 
Owen, Eilek, & Kassim 2000). For Sgr A* , we use the values M = 2.6 x 10^ Mq (e.g., 
Schodel et al. 2002) and M = 8.8 x 10"^ M© yr'^ (e.g.. Yuan, Markoff, & Falcke 2002; 
Quataert 2003). Although the kinetic luminosity of the jet in Sgr A* is rather uncertain (see, 
e.g., Yuan 2000; Yuan et al. 2002), we will adopt the value quoted by Falcke & Biermann 
(1999), and therefore we set Ljet = 5 x lO^^ergs s~^ 

We study both shocked and shock-free solutions spanning the computational domain 
between the inner radius Tin = 2.001 and the outer radius rout = 5000, where Tin and Tout 
are the radii at which the boundary conditions are applied (see eqs. [101]). Six different 
accretion/shock scenarios are explored in detail, with the values for the various parameters 
e_, i, e+, Tci, rc3, r*, fcs, H^, Al*, it!*, and T* reported in Table 1. Models 1, 2, and 3 are 
associated with M87, while models 4, 5, and 6 are used to study Sgr A* . In our numerical 
examples, we utilize natural gravitational units {GM — c = 1 and =2), except as noted. 
Based on the observational values for M and Ljet associated with the two sources, we can 
use equations (45) and (104) to conclude that Ae = —0.007 for M87 and Ae = —0.01 for 
Sgr A* . These results are consistent with the values for e_ and e+ reported in Table 1, and 
therefore Lghock = -^jet as required (sec cq. [104]). 

Next we use the energy conservation condition L^sc = Ljet (eq- [106]) to determine the 
value of the diffusion constant kq (eq. [70]) for a shocked disk. In Figure 3, we plot Lesc/Ljet, 
Fesc, and A^esc/-^o as functions of for the M87 and Sgr A* parameters, where Fesc = 
Eesc/ {fnpC^) is the mean Lorentz factor of the escaping particles. The treatment of energy 
conservation in our disk/shock model is self-consistent when the condition Lesc/Ljet = 1 is 
satisfied, which corresponds to specific values of kq as indicated in Figures 3a and 3d. We 
find that two kq roots exist for models 3 and 6, one root is possible for models 2 and 5, and 
no roots exist for models 1 and 4. Hence the values of e_ associated with models 2 and 5 
represent the maximum possible values for e_ that yield self-consistent solutions based on the 
M87 and Sgr A* data, respectively. For illustrative purposes, we shall focus on the details of 
the disk structure and particle transport obtained in models 2 and 5. 
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Fig. 3. — Quantities Lcsc/-^jct, Tcsc, and Ai"csc/^o for a shocked disk plotted as functions of 
kq. Panels (a), (6), and (c) correspond to the M87 parameters (models 1, 2 and 3), and 
panels (d), (e), and (/) correspond to the Sgr A* parameters (models 4, 5 and 6). The model 
numbers are indicated for each curve. See the discussion in the text. 



7.3. Disk Structure and Pcirticle Transport 



In order to illustrate the importance of the shock for the acceleration of high-energy 
particles, we shall examine the structure of the accretion disk either with and without a 
shock based on the values of the upstream parameters e_ and i utilized in models 2 and 5 
(see Table 1). In Figures 4a and 4& we plot the inflow speed v{r) and the adiabatic sound 
speed a(r) for the shocked and smooth (shock-free) solutions associated with models 2 and 5, 
respectively. Since we are working within the isothermal shock scenario, the sound speed a 
is continuous at the shock location. In Figures 4c and Ad we plot the specific internal energy 
U/p = {'-f — l)~^kT/mp for the background (thermal) gas along with the specific gravitational 
potential (binding) energy GM/{r — r^) as functions of radius for 7 = 1.5. These results 
demonstrate that the gas is marginally bound in the absence of a shock, and strongly bound 
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Fig. 4. — The velocity v{r) {solid lines) and sound speed a(r) [2/(7 + 1)]-'^/^ {dashed lines) are 
plotted in units of c for (a) model 2 and (6) model 5. The curves cross at the critical points. 
Also plotted are the specific internal energy U/p {solid lines) and the specific gravitational 
potential energy GM/{r — r^) {dashed lines), both in units of c^, for (c) model 2 and {d) 
model 5. The shocked and shock-free solutions are denoted by the thin and heavy lines, 
respectively 

when a shock is present. The increased binding of the thermal gas in the disk results from the 
escape of energy in the outflow, which reduces the sound speed compared with the shock-free 
case. The enhanced cooling allows the accretion to proceed, thereby removing one of the 
major objections to the original ADAF scenario (Narayan & Yi 1994, 1995). We emphasize 
that these new results represent the first fully self-consistent calculations of the structure of 
an ADAF disk coupled with a shock-driven outflow, hence extending the heuristic work of 
Blandford & Begelman (1999) and Becker, Subramanian, & Kazanas (2001). 

Next we study the solutions obtained for the relativistic particle number and energy 
density distributions in the disk based on the flow structures associated with models 2 and 
5. The related transport parameters are listed in Table 2. In Figures 5 and 6 we plot the 



-35- 



global number and energy density distributions obtained in a shocked disk using the model 2 
and 5 parameters, respectively. We also include the corresponding results obtained in a 
shock-free (smooth) disk for the same values of the upstream energy transport rate e_ and 
the specific angular momentum I. In each case the densities decrease monotonically with 
increasing radius. The increase near the horizon is a consequence of advection, while the 
decline as r — > oo reflects the fact that the particles injected at the shock have a very small 
chance of diffusing to large distances from the black hole. Note that the shocked disk has 
a lower value for the number density rij. at all radii as a consequence of particle escape. 
However, the shocked disk also displays a higher value for the energy density Ur, which 
reflects the central role of shock in accelerating the relativistic test particles. 

The kinks in the energy and number density distributions at the shock radius r — 
indicated in Figures 5 and 6 reflect the derivative jump conditions given by equations (75) 
and (95). The values for the ratios Nj/Nu and Nesc/No reported in Table 2 indicate that 
most of the injected particles are advected into the black hole, with ~ 20% escaping to form 
the outflow (sec Figs. 3c and 3/). In order to validate the accuracy of the numerical solutions 
for nr{r) and Ur{r), we also compare the proflles obtained with the asymptotic relations 
developed in § 6. We demonstrate in Figure 7 (model 2) and Figure 8 (model 5) that the 
solutions for both n,.(r) and Ur{r) agree closely with the asymptotic expressions given by 
equations (89) and (98) for small radii and by equations (86) and (96) for large radii. Note 
that the values reported by Le & Becker (2004) for n* and C/* were expressed in incorrect 
units and are given correctly in our Table 2. 



7.4. Jet Formation in M87 and Sgr A* 

The mean energy of the relativistic particles in the disk is given by (cf. eq. [107]) 

(E)^^, (109) 
nr{r) 

so that (E) = Eesc at r = r*. In Figure 9 we plot the mean energy as a function of radius 
in shocked and shock-free disks based on the parameters used for models 2 and 5. The 
results demonstrate that when a shock is present in the flow, the relativistic particle energy 
is boosted by a factor of ~ 5 — 6 at the shock location. By contrast, we flnd that in the 
shock-free models with the same values for e_, i, and Kq, the energy is boosted by a factor 
of only ~ 1.4 — 1.5. This clearly demonstrates the essential role of the shock in efficiently 
accelerating particles up to very high energies, far above the energy required to escape from 
the disk. Note that close to the event horizon, the mean energy of the relativistic particles is 
further enhanced by the strong compression of the accretion flow, as indicated by the sharp 
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Fig. 5. — Global solutions for the relativistic number density (a) and the relativistic energy 
density (b) computed using the model 2 parameters. The solid and dashed curves correspond 
to disks with and without shocks, respectively. Note that the number density is higher in 
the smooth (shock-free) disk due to the absence of particle escape. Conversely, the energy 
density is higher in the shocked disk due to the enhanced particle acceleration occurring at 
the shock. 

increase in (E) as r — > Tg . 

The material in the outflow is initially ejected from the disk in the vicinity of the shock 
as a hot plasma which cools as it expands, with its outward acceleration powered by the 
pressure gradient in the surrounding plasma. Based on our results for models 2 and 5, we 
find that the shock/jet locations are given by r* ~ 22 and r* ~ 16 for M87 and Sgr A* , 
respectively. The terminal (asymptotic) Lorentz factor of the jet, F^, can be estimated by 
writing 

Too = Tesc = ^ , (110) 

which is based on the assumption that the jet starts off "slow" and "hot" and subsequently 
expands to become "fast" and "cold." Adopting the Fesc values fisted in Table 2 for M87 
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Fig. 6. — Same as Figure 5 but for model 5. 

and Sgr A* , we obtain Poo = 7.92 (see Fig. 36) and Fo 



7.26 (see Fig. 3e), respectively. 



We can now compare our model predictions for the shock/jet location and the asymp- 
totic Lorentz factor with the observations of M87 and Sgr A* . According to Biretta, Junor, 
& Livio (2002), the M87 jet forms in a region no larger than ~ 30 gravitational radii from 
the black hole, which agrees rather well with our predicted shock/jet location ~ 22 for 
this source. Turning now to the asymptotic (terminal) Lorentz factor, we note that Biretta, 
Sparks, & Macchetto (1999) estimated > 6 for the M87 jet, which is comparable to the 
result Foo = 7.92 obtained using our model. In the case of Sgr A* , our model indicates 
that the shock forms at r* ~ 16 which is fairly close to the value suggested by Yuan (2000). 
However, future observational work will be needed to test our prediction for the asymp- 
totic Lorentz factor of Sgr A* , since no reliable observational estimate for that quantity is 
currently available. 
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Fig. 7. — Plots of the numerical solutions for ^^(r) and Ur{r) {solid lines) computed using 
the model 2 parameters in a shocked disk are compared with the asymptotic expressions 
{dashed lines) close to the event horizon (a and b) and at large radii (c and d). See the 
discussion in the text. 



7.5. Radiative Losses from the Jet 



It is still unclear whether the outflows observed to emanate from many radio-loud sys- 
tems containing black holes are composed of an electron- proton plasma or electron-positron 
pairs, or a mixture of both. Whichever is the case, the particles must maintain sufficient 
energy during their journey from the nucleus in order to power the observed radio emission, 
unless some form of reacceleration takes place along the way, due to shocks propagating along 
the jet (e.g., Atoyan & Dermer 2004a). Proton-electron outflows, such as those studied here, 
have a distinct advantage in this regard since most of the kinetic power is carried by the ions, 
which do not radiate much and are not strongly coupled to the electrons under the typical 
conditions in a jet (e.g., Felten 1968; Felten, Arp, & Lynds 1970; Anyakoha et al. 1987; 
Aharonian 2002). We therefore suggest that if the observed outflows are proton-driven, then 
they may be powered directly by the shock acceleration mechanism operating in the disk, 
with no requirement for additional in situ reacceleration in the jet. In this section we confirm 
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Fig. 8. — Same as Figure 7 but for model 5. 



this conjecture by considering the energy losses experienced by the protons in the outflow. 
The ions in the jet lose energy via two distinct channels, namely (1) direct radiative losses 
due to the production of synchrotron and inverse-Compton emission, and (2) indirect radia- 
tive losses via Coulomb coupling with the electrons. We will evaluate these two possibilities 
by estimating the corresponding coohng timescales for the outflows in M87 and Sgr A* and 
comparing the results with the jet propagation timescales for these sources. 

The energy loss rate due to the production of synchrotron and inverse-Compton emission 
by the relativistic protons escaping from the disk with mean energy Fesc^ipC^ is given by (see 
eqs. [7.17] and [7.18] from Rybicki & Lightman 1979) 

where Uph and Ub = B^/{8'k) denote the energy densities of the soft radiation and the 
magnetic held with strength B, respectively. The associated energy loss timescale is therefore 

t.^ - ^.72^! = ^ ( (Ub + t'pj-' . (112) 
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Fig. 9. — Mean energy of the relativistic particles in the disk, {E) = Ur{r) /n.f.{r) (eq. [109]), 
plotted in units of the injection energy £"0 as a function of radius for (a) model 2 and (6) 
model 5. Results are indicated for both shocked {solid lines) and shock-free {dashed lines) 
disk structures. 

In our apphcation to M87, we take B ~ 0.1 G based on estimates from Biretta, Stern, 
& Harris (1991), and we set Fesc ~ 8 (see Table 2, model 2). Assuming equipartition 
between the magnetic field and the soft radiation, this yields for the radiative cooling time 
^rad ~ lO^^yrs, which suggests that the protons can easily maintain their energy for many 
millions of parsecs without being seriously effected by synchrotron or inverse-Compton losses, 
as expected. For Sgr A* , we assume equipartition with i? ~ 10 G (Atoyan & Dermer 2004b) 
and Fesc ~ 7 (see Table 2, model 5). The radiative cooling time for the escaping protons 
is therefore t^ad ~ 10^ yrs. Hence synchrotron and inverse-Gompton losses have virtually no 
effect on the energy of the protons in the Sgr A* jet. 

In addition to synchrotron and inverse-Compton radiation, the protons in the jet will 
also lose energy due to Goulomb coupling with the thermal electrons, which radiate much 
more efficiently than the protons. The energy loss rate for this process can be estimated 
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using equation (4.16) from Mannheim & Schlickeiser (1994), which yields 

= 30 nga^meC^ , (113) 



T 

coul 



where rif. represents the electron number density in the jet. The associated loss timescale for 
a proton escaping from the disk with mean energy Fesc^^pC^ is 

^ Tcsc^pC rQsc (114) 

{dE / dt)\^^^^ SQuea^cme' 

The electron number density Ue decreases rapidly as the jet expands from the disk into the 
external medium. Hence the most conservative estimate (based on the strongest Coulomb 
coupling) is obtained by adopting conditions at the base of the jet, where rif. has its maximum 
value. To estimate the electron number density at the base of the outflow, we begin by 
calculating the rate at which protons escape from the disk at the shock location. By using 
equation (B8) to eliminate Aq in equation (76), we find that the proton escape rate is given 

N„,= '"-'y'\ (115) 

where r^,, n*, H^,, and Amag denote the radius, the proton number density, the vertical half- 
thickness, and the magnetic mean free path inside the disk at the shock location, respectively. 
The shock is expected to have a width comparable to Amag, and therefore the sum of the 
upper and lower face areas of the shock annulus is equal to 47rr*Ainag- We also note that the 
flux of the relativistic protons escaping from the disk into the outflow is given by cUp, where 
Up is the proton number density at the base of the jet. Combining these relations, we can 
write the proton escape rate in terms of Up using 

-^esc 47rrHcAmag 

CUp . (116) 

By equating the two expressions for A^esc given by equations (115) and (116), we find 
that Up is related to n* via 

!^^^ < 1. (117) 
if* 

Since the electron-proton jet must be charge neutral, the electron number density at the 
base of the jet, rie, is equal to the proton number density Up, and therefore we obtain 

rze = ^ n* . (118) 

Using the relation AmagZ-f^* — ^o^^ ^Q- [^8]) along with the results for Aq and reported 
in Table 1 for M87 (model 2), we obtain = 0.11 n* = 2.2 x lO^cm^^ Setting T^^ ~ 
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8, we find that equation (114) yields for the electron-proton Coulomb coupling timescale 
tcoui ~ 3.5 X lO^yrs. Note that this is an extremely conservative estimate since it is based 
on conditions at the bottom of the jet, and therefore it suggests that Coulomb coupling 
between the protons and the electrons is insufficient to seriously degrade the energy of the 
accelerated ions escaping from the disk as they propagate out to the radio lobes via the jet. 
For Sgr A* , we use the model 5 data in Table 2 to obtain rie = 0.13 = 5.6 x lO^cm"^. 
Setting Fesc ~ 7 yields for the Coulomb coupling timescale tcoui ~ 1.2 x 10^ yrs, which implies 
that the length of the jet can be as large several thousand parsecs before much energy is 
drained from the protons, assuming the material in the jet travels at half the speed of light. 
We emphasize that these numerical estimates of the importance of radiative and Coulomb 
losses experienced by the relativistic protons are based on the "worst-case" assumption that 
the conditions at the base of the outflow prevail throughout the jet. In reality, the jet density 
will drop rapidly as the gas expands, and therefore the true values for the proton energy loss 
timescales will be much larger than the results obtained here. This strongly suggests that 
shock acceleration of the protons in the disk, as investigated here, is sufficient to power the 
observed outflows without requiring any reacceleration in the jets. 



7.6. Radiative Losses from the Disk 

In the ADAF scenario that we have focused on, radiative losses from the disk are 
ignored. The self-consistency of this approximation can be evaluated by computing the 
free-free emissivity due to the thermal gas in the disk. The total X-ray luminosity can be 
estimated by integrating equation (5.15b) from Rybicki & Lightman (1979) over the disk 
volume to obtain for pure, fully-ionized hydrogen 

/■oo 

L,ad= / lA X 10-^^ T^/^ m;UV , (119) 

where dV = ATirHdr represents the differential (cylindrical) volume element, and Te denotes 
the electron temperature. We can obtain an upper limit on the X-ray luminosity by assuming 
that the electron temperature is equal to the ion temperature. Based on the detailed disk 
structures associated with models 2 and 5, we find that Lj.ad/-^jct ~ 10~^ and Lrad/-^jct ~ 
10~^, respectively. However, in an actual ADAF disk, the X-ray luminosity will of course 
be substantially smaller than these values because the electron temperature is roughly three 
orders of magnitude lower than the ion temperature. Hence our neglect of radiative losses is 
completely justified, as expected for ADAF disks. 
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8. CONCLUSION 

In this paper we have demonstrated that particle acceleration at a standing, isothermal 
shock in an ADAF accretion disk can energize the relativistic protons that power the jets 
emanating from radio-loud sources containing black holes. The work presented here rep- 
resents a new type of synthesis that combines the standard model for a transonic ADAF 
flow with a self-consistent treatment of the relativistic particle transport occurring in the 
disk. The energy lost from the background (thermal) gas at the isothermal shock loca- 
tion results in the acceleration of a small fraction of the background particles to relativistic 
energies. One of the major advantages of our coupled, global model is that it provides a 
single, coherent explanation for the disk structure and the formation of the outflow based 
on the well-understood concept of first-order Fermi acceleration in shock waves. The the- 
ory employs an exact mathematical approach in order to solve simultaneously the combined 
hydrodynamical and particle transport equations. 

The analysis presented here closely parallels the early studies of cosmic-ray shock accel- 
eration. As in those first investigations (e.g., Blandford & Ostriker 1978), we have employed 
an idealized model in which the pressure of the accelerated particles is assumed to be negligi- 
ble compared with that of the thermal background gas (the "test particle" approximation) . 
In order to check the self-consistency of this assumption, we have confirmed that the total 
pressure is dominated by the pressure of the background (thermal) gas throughout most 
of the disk. However, in the vicinity of the shock the two pressures can become compa- 
rable and this suggests that the dynamical results will change slightly if the test particle 
approximation is relaxed. We plan to consider this question in future work by developing 
a "two-fiuid" version of our model that includes the particle pressure, in analogy with the 
"cosmic-ray modified shock" scenario for cosmic-ray acceleration (Becker & Kazanas 2001; 
Drury & Volk 1981). 

We have presented detailed results that confirm that the general properties of the jets 
observed in M87 and Sgr A* can be understood within the context of our disk/shock/outfiow 
model. In particular, our results indicate that the shock acceleration mechanism can produce 
relativistic outflows with terminal Lorentz factors and total powers comparable to those 
observed in M87 and Sgr A* . However, in principle even higher efficiencies can be achieved 
by varying the upstream energy transport rate e_ which is the fundamental free parameter 
in our model. The buildup of the particle pressure in such high-efficiency situations would 
require relaxation of the test-particle approximation, as discussed above. In this paper 
we have focused on inviscid disks, which are the simplest to analyze. While the inviscid 
model provides useful insight into the importance of shock acceleration in ADAF disks, this 
restriction clearly must be lifted in the future, since viscosity plays a key role in determining 
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the structure of an actual accretion disk. We are currently developing a self-consistent 
viscous disk model in order to explore shock formation and particle acceleration in a more 
rigorous context. However, we do not expect the presence of viscosity to alter any of the 
basic conclusions reached in this paper because significant particle acceleration will occur 
regardless of the viscosity, provided a shock is present. The existence of shocks in viscous 
disks is a controversial issue, but several studies suggest that shock formation is possible 
provided the viscosity is relatively low. In the absence of a consensus regarding the possible 
presence of shocks in accretion disks, we believe that it is important to study models with 
shocks in order to develop theoretical predictions that can be tested observationally. 

The shock acceleration mechanism analyzed in this paper is effective only in rather 
tenuous, hot disks, and therefore we conclude that our model may help to explain the 
observational fact that the brightest X-ray AGNs do not possess strong outflows, whereas 
the sources with low X-ray luminosities but high levels of radio emission do. We suggest 
that the gas in the luminous X-ray sources is too dense to allow efficient Fermi acceleration 
of a relativistic particle population, and therefore in these systems, the gas simply heats 
as it crosses the shock. Conversely, in the tenuous ADAF accretion flows studied here, the 
relativistic particles are able to avoid thermalization due to the long collisional mean free 
path, resulting in the development of a significant nonthermal component in the particle 
distribution which powers the jets and produces the strong radio emission. We therefore 
conclude that the coupled, self-consistent theory for the disk structure and the particle 
acceleration investigated here is capable of powering the outflows observed in many radio- 
loud systems containing black holes. 

The authors are grateful to Dr. Lev Titarchuk for providing a number of useful com- 
ments on the manuscript, and also to the anonymous referee for several insightful suggestions 
that signiflcantly improved the paper. 
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APPENDICES 

A. Treatment of the Vertical Structure 

In principle, the pressure P, density p, diffusion coefficient k, Green's function f^, and 
velocity components Vr and in the disk all display significant variations in the vertical (z) 
direction. Following Abramowicz & Chakrabarti (1990), we will use the first five quantities to 
represent vertical averages over the disk structure at radius r. However, the vertical variation 
of the velocity component must be treated differently. Here, we assume for simplicity that 
the vertical expansion is homologous, and therefore the vertical velocity variation is given by 

v,{r,z)^B{r)z . (Al) 

It follows that the vertical velocity at the surface of the disk, z — H{r), can be written as 



Vzir,z) 



z=H 



B{r)H{r) . 



(A2) 



In a steady-state situation, we can also express the vertical velocity at the disk surface using 

dH 



Vz{r,z) „ 

z=H dr 

By combining the two previous expressions, we find that the function B{r) is given by 

dlnH 

B{r) = Vr 



dr 



(A3) 



(A4) 



This result will prove useful when we vertically integrate the transport equation. Note that 
in terms of S(r), we can write the divergence of the fiow velocity v in cylindrical coordinates 
as 

- Id , . dvz Id ^ . ^^^^ 



V • -y = {rvr) + ^ = {rvr) + B{r) , 
r or oz r or 

where we have assumed azimuthal symmetry. Application of equation (A4) now yields 

^■'^^WrWr^'^'''^ ■ ^^^^ 



The steady-state transport equation expressed in cyhndrical coordinates is (see eq. [61]) 



dr 



dz 



Id, . dvz 
r or dz 



dE r dr \ dr 



+ 



NoS{E-Eo) S{r-r,) 
(AnEoYr^H, 



Aocd{r - r*) . 



(A7) 
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Operating on equation (A7) with dz and applying equation (Al) yields, after partially 
integrating the term containing Vz on the left-hand side, 



d 



1 d 
r dr 



[rvr) + B 



HE 



dE r dr \ 



i tHk 



dr 



+ 



N^5{E-Eo)5{r-r,) 



AocH,5{r-r,)f^ , (A8) 



(47r£;o)2r. 

where the symbols /q, Vr, and k now refer to vertically averaged quantities. Using equa- 
tions (A4), (A5), and (A6), we can rewrite the vertically integrated transport equation as 



Hv, 



= —TT {fHv,) E^ + - — 
or 3r or oE r or 

NoS{E-Eo)S{r-r,) 



rHn 



dfa 
dr 



+ 



i^nEoyr, 



AocH^S{r - r^) . 



(A9) 



B. Derivation of the Escape Pcirameter 

The dimensionless parameter Aq appearing in equation (58) determines the rate of 
particle escape through the surface of the disk due to random walks occurring near the shock 
location. Since the particles are accelerated as a consequence of collisions with magnetic 
waves, we will assume that the thickness of the shock is comparable to the magnetic mean 
free path, Amag- In order to estimate Aq, we model the escape of the particles from the disk 
using the analogy of "leakage" through an opening in a cylindrical pipe with radius equal to 
the half-thickness of the disk at the shock location, if*. The length of the open section of 
the pipe is set equal to the shock thickness Amag- The particle number density in the open 
section is governed by the equation 

drij- Tij- /T->i \ 

Vx—r- — , Bl 

dr t 

where v^, rir, and tcsc represent the flow velocity, the relativistic particle number density, and 
the average time for the particles to escape through the open walls of the pipe via diffusion. 
Upon integration, the solution to equation (Bl) is given by 

nr{x) = no exp ( ^ J , (B2) 

where no is the incident number density as the flow encounters the opening in the pipe, at 
a; = 0. We can approximate the solution for nr{x) by performing a Taylor expansion around 
X = 0, which yields 

n,(x) ft; no ( 1 - ^— ) . (B3) 
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The fraction of particles that escape from the pipe can therefore be estimated by setting 
X — '^mag to obtain 

/esc = l-- = ^. (B4) 

In order to make contact with the disk apphcation, we note that according to equations (75) 
and (76), the fraction of particles that escape as the gas crosses the isothermal shock is given 

/esc = ^ — , (B5) 

where v^, = (f+ + ?;_)/2 is the mean velocity at the shock, and we have assumed that 
advection dominates over diffusion. Eliminating /esc between equations (B4) and (B5), and 
setting Vx = v^, we find that 

Ao^^. (B6) 

^ ^esc 

Within the context of our one-dimensional model for the particle transport in the disk, 
the mean escape time tesc is related to Amag and the disk half-thickness at the shock if* via 

H 

tesc = = , * , (B7) 

where v^m — cAmag/ denotes the vertical diffusion velocity of the protons in the tangled 
magnetic field near the shock, which is valid provided if/ Amag > 1- Eliminating tesc between 
equations (B6) and (B7) then yields 

^0= (^)' < 1. (B8) 

The diffusion coefficient at the shock is related to the magnetic mean free path by the 
standard expression (e.g., Reif 1965) 

(B9) 

and therefore equation (B8) can be rewritten as 

^„^(^)\ (BIO) 

where = (/t_ + /t+)/2 denotes the average of the upstream and downstream values of k 
on either side of the shock. 
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